Introducion
The Diel.Niche package is used to evaluate how animals use the three
fundamental periods of light: twilight (dawn & dusk), daytime, and
nighttime. First, install and load the library
#Install package from simplified repo.
#Make sure to add your auth_token
#devtools::install_github("diel-project/Diel-Niche-Modeling"
# ,ref="Diel.Niche.Simplify"
# ,auth_token = "")
#Call the library
library(Diel.Niche)
Diel.Niche specifies many diel hypotheses, which together are defined
into five hypothesis sets. The code names for the hypotheses and sets
are,
what.hyp()
#> Traditional : D N CR C
#>
#> General : D N CR C2 D.CR D.N CR.N
#>
#> Threshold : D.th N.th CR.th C.th
#>
#> Maximizing: D.max N.max CR.max
#>
#> Variation: D.var N.var CR.var C.var
#>
A short description of each hypothesis is available by calling,
what.hyp("D")
#> [1] "Traditional Diurnal"
To list the hypotheses of a given set,
hyp.sets("Traditional")
#> [1] "D" "N" "CR" "C"
hyp.sets("Threshold")
#> [1] "D.th" "N.th" "CR.th" "C.th"
The best way to understand each hypothesis set is by plotting all
hypotheses.
Maximizing
This hypothesis set is most useful if you are only interested in
evaluating which diel modality (Diurnal, Nocturnal, Crepuscular) is used
the most. There is no hypothesis associated to cathemerality, which is a
combination of multiple modalities.
diel.plot(hyp=hyp.sets("Maximizing"))
Traditional
This hypothesis set uses the traditional combinations of the diel
modalities of diurnal, nocturnal, crepuscular, and cathemeral.
diel.plot(hyp=hyp.sets("Traditional"))
If left unspecified, default values are used to define the boundaries
among these different hypotheses. The default value for being diurnal,
nocturnal, or crepuscular is \(\ge\)
0.8 probability. Cathemerality is defined as no modality has more than
0.8 probability and two or more modalities have\(\ge\) 0.1 probability.
General
This hypothesis set includes the Traditional definitions of diurnal,
nocturnal, and crepuscular, as above, but cathemerality is defined more
specifically in terms of whether two or three modalities have a high
probabilty.
diel.plot(hyp=hyp.sets("General"))
Threshold
This hypothesis set is based on defining each modality in terms of
lower thresholds. If the thresholds are not provided, default values
ared used.
diel.plot(hyp=hyp.sets("Threshold"))
Variation
This hypothesis set is based on defining the ranges (lower, uppper)
of possible probability values.
diel.plot(hyp=hyp.sets("Variation"))
We can also view specific hypotheses, as
diel.plot("D",more.points = TRUE)
Setup hypothesis set and fit models
You can setup your own model set or use the pre-defined model
sets.
hyp.set = c("D","CR") #manual
hyp.set = hyp.sets("Traditional") #pre-defined
To fit the data to the hypothesis set,
out = diel.fit(y,hyp.set)
#> Data checks Complete.
#> Calculating Bayes Factors...
#> The most supported model is:
#> CR
#Call the model probabilities for each hypothesis in the set
out$bf.table
#> Prior Posterior
#> D 0.25 0.0000000000
#> N 0.25 0.0000000000
#> CR 0.25 0.9992042477
#> C 0.25 0.0007957523
To examine the available outputs from the fitted model object,
attributes(out)
#> $names
#> [1] "bf.table" "post.samp" "ms.model"
#> [4] "ppc" "ms.ppc" "post.samp.ms.model"
#> [7] "y" "y.vec" "gelm.diag"
#> [10] "ms.gelm.diag" "bf.list"
?diel.fit
#> starting httpd help server ... done
The function ‘diel.fit’ defaults to providing the model probabilities
for each hypothesis set, but not the posterior samples of the parameters
for each hypothesis. We can change this as,
out = diel.fit(y,hyp.set,n.chains=2,post.fit = TRUE)
#> Data checks Complete.
#> Calculating Bayes Factors...
#> Posterior Sampling...
#> The most supported model is:
#> CR
To look at convergence criteria for each hypothesis. Note that
convergence can be very poor for models that do not fit the data
well.
out$gelm.diag
#> $D
#> Potential scale reduction factors:
#>
#> Point est. Upper C.I.
#> p1_1 0.999 1
#> p1_2 1.001 1
#>
#>
#> $N
#> Potential scale reduction factors:
#>
#> Point est. Upper C.I.
#> p1_1 NaN NaN
#> p1_2 1 1.01
#>
#>
#> $CR
#> Potential scale reduction factors:
#>
#> Point est. Upper C.I.
#> p1_1 1.003 1.02
#> p1_2 0.999 1.00
#>
#>
#> $C
#> Potential scale reduction factors:
#>
#> Point est. Upper C.I.
#> p1_1 1 1.00
#> p1_2 1 1.02
We can plot the posterior samples from the most supported model to
check convergence/mixing as,
plot(coda::as.mcmc(out$post.samp.ms.model))
The posterior samples for all hypotheses are available in a list.
names(out$post.samp)
#> [1] "D" "N" "CR" "C"
#For each of these list is a list of chains
length(out$post.samp[[1]])
#> [1] 2
#Here are the means of the posterior samples of all hypotheses for chain 1
lapply(out$post.samp,FUN=function(x){colMeans(x[[1]])})
#> $D
#> p_crep_1 p_day_1 p_night_1
#> 0.17646753 0.80201824 0.02151423
#>
#> $N
#> p_crep_1 p_day_1 p_night_1
#> 0.0000000 0.1061735 0.8938265
#>
#> $CR
#> p_crep_1 p_day_1 p_night_1
#> 0.87862454 0.01813177 0.10324369
#>
#> $C
#> p_crep_1 p_day_1 p_night_1
#> 0.78413959 0.03457643 0.18128398
Plotting
Using the packages bayesplot and ggplot2, we can examine our
posterior distributions along with the true probabilities values,
library(ggplot2)
library(bayesplot)
#> This is bayesplot version 1.9.0
#> - Online documentation and vignettes at mc-stan.org/bayesplot
#> - bayesplot theme set to bayesplot::theme_default()
#> * Does _not_ affect other ggplot2 plots
#> * See ?bayesplot_theme_set for details on theme setting
posteriors=coda::as.mcmc(out$post.samp.ms.model)
plot_title <- ggtitle("Posterior distributions",
"with medians and 80% intervals")
mcmc_areas(posteriors, prob = 0.8) + plot_title+
geom_vline(xintercept=p, linetype="dashed",
color = c("red","green","purple"), size=1)
Another way to examine our hypothesis is to plot the theoretical
niche space for the hypothesis along with the posterior samples from the
most supported model.
diel.plot(hyp=out$ms.model,more.points = TRUE,
posteriors=out$post.samp.ms.model)
Or, we can plot the whole hypothesis set
diel.plot(hyp=hyp.set,
posteriors=out$post.samp.ms.model)
User-specified Hypotheses
You do not need to rely on pre-set definitions of diel hypotheses. If
you are interested in categorizing a diel modality based on 0.9
probability or higher, you can adjust the hypotheses using the function
‘diel.ineq()’ than pass this object to the plot function or model
fitting function.
diel.setup=diel.ineq(xi.min.dom = c(0.9))
diel.plot(hyp=hyp.sets("Traditional"),diel.setup=diel.setup)
If you want to separate the hypotheses and reduce the parameter space
of cathmerality,
diel.setup=diel.ineq(xi.min.dom = 0.7,separation=0.1)
diel.plot(hyp=hyp.sets("Traditional"),diel.setup=diel.setup)
When considering the General hypothesis set, you can adjust the lower
probability value to define crepuscular, diurnal, and nocturnal, as well
as adjust the definitions of Cathemeral General vs Diurnal-Nocturnal,
Crepuscular-Nocturnal, and Diurnal-Crepuscular.
diel.setup=diel.ineq(xi.min.dom = c(0.9,0.05))
diel.plot(hyp=hyp.sets("General"),diel.setup=diel.setup)
This new hypothesis set can be used to simulate data and fit the same
models by passing the object diel.setup as an attribute of the functions
sim.diel, diel.fit, and diel.plot.
y=sim.diel(hyp="D.CR",n.sample = 200,diel.setup=diel.setup)$y
out = diel.fit(y,hyp.set=hyp.sets("General"),n.chains=2,post.fit = TRUE,diel.setup=diel.setup)
#> Data checks Complete.
#> Calculating Bayes Factors...
#> Posterior Sampling...
#> The most supported model is:
#> D.CR
diel.plot(hyp=hyp.sets("General"),posteriors = out$post.samp.ms.model,diel.setup=diel.setup)